Interdiffusion of Solvent into Glassy Polymer Films: A Molecular 

Dynamics Study 

Mesfin Tsige* and Gary S. Gresv 
Sandia National Laboratories, Albuquerque, NM 87185 
(Dated: February 2, 2008) 

Abstract 

Large scale molecular dynamics and grand canonical Monte Carlo simulation techniques are 
used to study the behavior of the interdiffusion of a solvent into an entangled polymer matrix 
as the state of the polymer changes from a melt to a glass. The weight gain by the polymer 
increases with time t as i 1 / 2 in agreement with Fickian diffusion for all cases studied, although 
the diffusivity is found to be strongly concentration dependent especially as one approaches the 
glass transition temperature of the polymer. The diffusivity as a function of solvent concentration 
determined using the one-dimensional Fick's model of the diffusion equation is compared to the 
diffusivity calculated using the Darken equation from simulations of equilibrated solvent-polymer 
solutions. The diffusivity calculated using these two different approaches are in good agreement. 
The behavior of the diffusivity strongly depends on the state of the polymer and is related to the 
shape of the solvent concentration profile. 
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I. INTRODUCTION 



Understanding the interdiffusion of solvent into a polymer is crucial for a variety of ap- 
plications, such as food storage, controlled drug release and membrane separations.— i 2 i 3 i 4 i 5 
The mechanisms controlling the interdiffusion process are reasonably well understoo d 6 ^ 8 ' 9 
but predicting accurately the nature of the diffusion has been a challenging problem. It 
is now widely accepted that the interdiffusion of solvent into a polymer depends on sol- 
vent concentration gradient within the system as well as the rate of polymer segmental 
relaxation . 10 i 11 i 12 i 13 i 14 Whether the polymer is a melt above its glass transition temperature 
T g or an amorphous solid below T g strongly affects the diffusion behavior.— 

In general, three categories of diffusion behavior of solvents into polymers have been 
distinguished. ^ 1 ^ 16 These are, Fickian or Case I, Case II or Class II, and anomalous dif- 
fusion: in which the rate of diffusion of solvent is much less than, much greater than, or 
comparable to the rate of polymer segmental relaxation, respectively. A simple descriptive 
way of quantifying these is based on the power law dependence of the mass uptake of the 
solvent by the polymer or the distance covered by the solvent as a function of time t (~ t n ). 
For Fickian diffusion n = 1/2, for Class II diffusion n = 1, and for anomalous diffusion 
1/2 < n < 1. Fickian diffusion usually applies for all solvent concentration when the solvent 
interdiffuses into a polymer melt, while for glassy polymers it usually applies only for low 
solvent concentration. Non-Fickian kinetics is expected when the viscoelastic properties of 
the system becomes the determining factor . 10 i n In addition to linear kinetics, Case II diffu- 
sion is characterized by a sharp concentration front that propagates at constant spee d 11 ^ 7 ! 18 
with a Fickian type precursor foo t 8 ' 15 ! 18 ' 19 preceding the front. 

When a solvent film is placed in contact with one surface of a polymer melt, the diffusion 
is one-dimensional and can often be described by Fick's one-dimensional diffusion equation 20 



where c is the solvent concentration in units of mass per unit volume and D(c) is the 
diffusivity. This equation assumes that the volume of the medium is not changed by the in- 
terdiffusion of the solvent. If D(c) is a function of c only, then the Boltzmann transformation 
of Eq. [T] gives 
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where f(D, c) is a function of D and c only. This equation reflects the square root time de- 
pendence of Fickian diffusion irrespective of the functional form of D(c). It can be integrated 
to yield the diffusion coefficient at concentration c— 



where 77 = z/t 1 ^ 2 . Thus from the scaled concentration profile one can directly obtain the 
diffusivity D(c). Note that only in a few special cases, like D(c) constant, can Eq. ^ be 
solved analytically.-^ For the case in which D(c) is a constant, c(r]) is an error function. 
However, it is important to remember that Fickian diffusion, i.e. uptake increasing as £ 1//2 , 
is true provided that the diffusivity depends only on c as shown above. 

The diffusivity can also be approximately obtained from the following Darken equation 
applied to solvent diffusing in an equilibrated polymer solution 



where D c (c) is the corrected diffusion constant and / is the fugacity of the solvent, both are 
defined in the next section. 

Molecular dynamics (MD) simulation technique is proven to be a useful tool for determin- 
ing the diffusion coefficients of penetrant molecules in polymers. This technique is specially 
important when detailed microscopic information of the mechanism of transport is required. 
Most of the previous studies have focused on the penetrant transport of small molecules in a 
polymer me it 3 5 i 21 i 22 i 23 i 24 i 25 i 26 i 27 i 28 i 29 i 30 i 31 With recent advances in parallel molecular dynam- 
ics algorithms and the increased computational power, progress has primarily occurred for 
studying the diffusion of large molecules (phenol molecules) in a polymeric matrix at atom- 
istic level . 31132 However, equivalent development is lacking for interdiffusion of solvent into 
polymer or polymer-polymer interdiffusion. In the previous study, which will be referred to 
hereafter as paper I,— we investigated the interdiffusion of a solvent into a homopolymer 
melt. The solvent concentration profile and weight gain by the polymer was measured as a 
function of time. The weight gain was found to scale as t 1//2 and the concentration profiles 
were found to fit very well assuming Ficks's second law with constant diffusivity. The study, 
however, focused only on homopolymers that are far above the glass transition temperature. 

In this paper we extend our previous study on interdiffusion of solvent into homopolymers 
that are close to the glass transition temperature of the homopolymer. Case II diffusion has 
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been inaccessible to computer simulation due to the extensive computational effort required. 
In this paper we study the conditions that may lead to Case II type diffusion behavior. The 
main purpose of the present study is two-fold: first, to understand how the dependence 
of the diffusivity D(c) on concentration changes as the state of the polymer changes and 
also its relation to the form of the concentration profile curve; second, to test the Darken 
approach under general conditions where D c (c) is not a constant. As in paper I, we are 
also interested in the relationship between the self-diffusion constant of the solvent and the 
corrected diffusion and the diffusivity. 

The outline of this paper is as follows. In Sec. II a brief review of the molecular dynamics 
simulation and the model used is given. In Sec. Ill the interdiffusion results for different 
polymer-polymer and solvent-polymer interaction parameters are presented and discussed. 
The diffusivity D(c) calculated from solvent concentration profiles and from the Darken 
equation are presented in Sec. IV. The self- and corrected diffusion constants as a function 
of solvent concentration are also presented and discussed. The main results of the present 
study are summarized in Sec. V. 

II. SIMULATION DETAILS 
A. System 

The basic model of the polymer- solvent system is the same as used in paper I. The 
polymer is treated as freely jointed bead-spring chain of length N monomers of mass m 
and the solvent is modeled as single monomer of mass m. All monomers of type a and (3 
interacts through the standard Lennard- Jones 6-12 potential 



where r is the distance between monomers and eu is a constant added so that the potential 
is continuous at r = r c . Here a = p stands for the polymer monomer and a = s for a solvent 
monomer. e ss = e and a define the units of energy and length, respectively. Here we take 
o = Cap and r c = 2.5a. For our model, the freezing temperature of the solvent is higher than 
the glass transition temperature of a long fully flexible polymer melt, T g = 0.5 — 0.6e//ce-~ 
Thus, temperature is not a good variable to change the state of the polymer without changing 
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the state of the solvent. Instead, we vary e pp from e pp = e (melt) to 2e (glassy). Berthelot 
rule e sp = y/e ss e pp is used for the cross term in some cases but we also study other cross 
terms. This is because (e pp , e sp ) = (2e, v2e) is found to be immiscible except in the dilute 
limit. As in paper I, for bonded monomers an additional anharmonic potential known as 
FENE potential with Rq = 1.5a and k = 30e is applied.^^ 

For the interdiffusion of solvent monomers into a polymeric matrix, the system consist of 
entangled polymer chains in a rectangular box which is periodic in x and y but not in z, the 
diffusion direction. This initial configuration was generated following the procedure given 
in paper I. The polymer consisted of 600 chains of length N = 500 monomers. The solvent 
consisted of 230,000 monomers. For the self- and corrected diffusion studies as a function 
of concentration, the system consist of an equilibrated polymer solvent mixture in a cubic 
box which is periodic in all directions. The polymer in the system consisted of M chains 
of length N = 500 monomers, where M = 100 for solvent concentration c < 0.45a -3 and 
M = 50 for c > 0.45a -3 . The mole fraction of solvent x s in the mixture was varied from 
0.01 (dilute case) to 0.75. A pure solvent system of 50,000 monomers was also simulated. 

In paper I we compared results from Langevin thermostat simulations which screens the 
hydrodynamic interactions with results from dissipative particle dynamics (DPD) thermo- 
stat simulations, which does not. The results from the two thermostats agreed when the 
dissipation from the thermostats become much smaller than from particle collisions. In the 
present study, to conserve hydrodynamic interactions, we use DPD thermostat through out 
the simulation. For details see paper I. The equations of motion were integrated with a 
velocity verlet algorithm with a time step of At = 0.012r for the interdiffusion study and 
At = 0.009r for the bulk equilibrium measurement of the self- and corrected diffusion con- 
stants, where r = m(a / 'e) 1 / 2 . All the simulations were run using the massively parallel code 
LAMMPS 37 at a temperature of T = e/fcs and pressure P ~ without tail corrections to 
be comparable to interdiffusion - same as paper I. 

B. Diffusion Coefficients 

We have calculated the self-and corrected diffusion coefficients and diffusivity D(c) of 
solvents in an equilibrated solvent polymer mixture as a function of solvent concentration, 
c. The self-diffusion constant D s (c) of the solvent in the polymer was calculated from the 
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slope of the solvent mean square displacement 

^>^ <[r( V° )12> <«) 

Here (...) denote an ensemble average and is obtained by averaging over all solvents and 
many initial time origins, and r(t) is the position of a solvent in the polymer at time t. 

The corrected diffusion coefficient D c (c) of the solvent in the polymer is calculated using 
the Einstein form equation 33 

D c (c) = N T x s x p lim ^-({[r cmjS (t) - r cm;S (0)] 

-[rcm, P (t) -r cm , P (0)]} 2 > (7) 

where x« and r cm ^(t) are mole fraction and center of mass of all monomers of species i at 
time t, respectively, and Nt = N s + N p is the total number of monomers. 

To determine the fugacity / the particle insertion method 3 * 8 is applied using the grand 
canonical MD code LADERA. 39 During the course of an equilibrium molecular dynamics 
simulation at the appropriate solvent concentration, the energy, E, of inserting a solvent par- 
ticle at random locations was sampled. The excess chemical potential energy fi e is computed 
using 

fi e = -k B T\n(exp(-E/k B T)}, (8) 

where k B is the Boltzmann constant, T is the temperature and (...) is an ensemble average. 
Then, the activity coefficient, 7, is computed via 7 = exp(fi e /k B T). The thermodynamic 
factor in eq. 0] can be expressed in terms of the activity coefficient 7 of the solvent as 

9 In/ _ 9 In 7 , . 

o In c a In c 

The thermodynamic factor goes to 1 as c — > 0. 

The computing time depends on the state of the polymer where much longer run are 
required as the effective temperature of the polymer melt is reduced towards its glass tran- 
sition temperature. At the lowest temperature studied and for a given solvent concentration 
c a run of about half a million MD time steps are required to calculate the self-diffusion 
constant D s (c) while a run of more than four million MD time steps is required for the 
corrected diffusion constant D c (c). The fugacity calculation at each solvent concentration 
requires more than four million MC insertion attempts. For the interdiffusion simulations, 
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one run was made for each set of parameters specified. Each system is run until the solvent 
reaches the lower substrate and on the average requires a run of about five million MD time 
steps. 

III. INTERDIFFUSION 

Interdiffusion studies of a solvent into an equilibrated polymer has been conducted for 
different cases of polymer-polymer and solvent-polymer interactions. The initial setup for 
the interdiffusion study is the same as Fig. 1 of paper I. The density profile of both polymer 
and solvent as a function of time for k B T /e pp = 1.0, 0.75 and 0.5 with Berthelot's rule for the 
cross-term e sp is shown in Fig. ^ The solvent diffuses into the polymer from the right side. 
As the state of the polymer changes from melt (/c^T/e^ = 1.0) to glass {k B T/e pp = 0.5) 
the solvent density profile changes to a sharp front. The solvent diffusion is Fickian in all 
the three cases as confirmed by the linearity of the weight gain by the polymer versus t 1 ^ 2 
curve shown in Fig. 121 This indicates that the precursor of the front for k B T /e pp = 0.5 
is Fickian. This is in agreement with recent experimental observations that characterize 
Case II diffusion by a sharp concentration front with a Fickian type precursor . 8 ! 15 ' 18 ' 19 For 
the front to move the solvent mobility should be much greater than the rate of polymer 
segmental relaxation. 10 But, for this case the front does not move in the time scale of our 
simulation. In fact, simulation of an equilibrated solvent-polymer solution, discussed in the 
next section, shows that only a small amount of solvent is soluble for this case suggesting 
that the front may not move at all. 

In order to facilitate the interdiffusion of solvent into the glassy polymer (k B T '/e pp = 0.5), 
the interaction between polymer and solvent is increased to e sp = 1.55e, 1.7e, and 2.0e. 
Similarly e sp is increased to 1.33e for k B T ' /e pp = 0.75. The corresponding density profiles of 
both polymer and solvent as a function of time is shown in Fig. El We clearly see that as 
e sp is increased the solubility is enhanced for both e pp = 1.33e and 2.0e and the density front 
observed for (e pp , e sp ) = (2.0e, \/2X)e) disappears. The e sp = 1.7e case shows a cross-over. 
There is no change in the diffusion process due to the change in e sp as the weight gain by 
the polymer system for all four cases increases as t 1//2 , see Fig. in agreement with Fickian 
diffusion. 
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FIG. 1: Solvent and polymer concentration profiles as a function of time for (a) ksT '/e pp = 1.0, 
plotted every 2400t, (b) ksT/epp = 0.75, plotted every 4800r, and (c) ksT/e pp = 0.5, plotted 
everyl2000r, with Berthelot's rule e sp = ^/e pp e ss for the cross-term. The solvent diffusing into the 
polymer from the right side. 

IV. DIFFUSION COEFFICIENTS 

The diffusivity of the solvent can be calculated from the solvent concentration profile using 
Eq. 01 It can be also calculated from simulation of equilibrated solvent-polymer solution 
using the Darken equation (Eq. |3J). In this section we compare diffusivity results from the 
two different approaches. 

A. Diffusivity from concentration profile 

Using the change of variable i] = zt~ x l 2 the solvent density profiles corresponding to 
different times superimpose as shown in Fig. This indicates that the solvent diffusivity 
is independent of position as expected for Fickian diffusion. It also means that the solvent 
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FIG. 2: Weight gain in terms of the number of solvent monomers diffused into the polymer at a 
given time t, for (e pp , e sp )= (e, e) (circles), (1.33e, \/L33e) (squares), (2.0e, y/2.0e) (down triangles), 
(2.0e, 1.55e) (stars), (2.0e, 1.7e) (diamond), and (2.0e, 2.0e) (up triangles), T = e/k B . 

diffusivity for absorption can be determined as a function of concentration by integrating 
Eq. El with respect to rj. 

To calculate D(c) analytically using Eq. El the average of the transformed solvent density 
profiles c(rj) of each case shown in Fig. 0]was fit to a polynomial function of at least order 5. 
This higher order polynomial was chosen to get an optimal fit to the concentration profiles. 
The data is integrated analytically up to the target concentration using the transformation 
J ° rjdc = f v o rj^-dr]. This procedure is repeated for different values of solvent concentration 
and the diffusivity D(c) calculated for the different cases is shown in Fig. E3 

In general, the behavior of D(c) strongly depends on the state of the polymer. The diffu- 
sivity is approximately a constant when the homopolymer is far above the glass transition, 
T g , that is for kBT/e pp = 1.0, and then becomes concentration dependent as the effective 
temperature of the polymer melt fcsT/e pp is reduced towards its glass transition tempera- 
ture. For (e pp , e sp ) = (1.33e, a/1.33c), (2.0e, v2\0e) and (2.0e, 1.7e) the calculated diffusivity 
for concentrations in the steep part of the concentration profile (not included in Fig. 03) re- 
sulted in large scatter of the data. This is because little variation of the slope in this region 
introduces large error in the diffusivity. 

For e pp = 1.33e, the diffusivity is approximately a constant for e sp = vL33e, but increases 
linearly with concentration for e sp = 1.33e. Note that in the dilute limit (c — > 0) the 
diffusivity is independent of e sp . For e pp = 2.0e, the diffusivity is independent of e sp within the 
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FIG. 3: Solvent p s and polymer p p concentration profiles as a function of time, plotted every 
2400r for (e pp , e sp ) = (a) (1.33e,1.33e), (b) (2.0e,1.55e), (c) (2.0e,1.7e), and (d) (2.0e,2.0e). 

error of the simulation. For this case the diffusivity can be approximated by an exponential 
function of the form D(c) = D exp(ac), where D and a are constants that may depend on 
the state of the polymer, and is represented by a line in Fig. |5^b). The empirical expression 
found from fitting the diffusivity curve of a given system was in turn used to solve Eq. Q 
and the calculated concentration profiles are shown in Fig. 0] as solid lines. In all cases, the 
calculated concentration profiles give an adequate description of the simulated concentration 
profiles in the region of interest. 

The behavior of D(c) is directly related to the form of the concentration profile 
curve. When the solvent concentration profile is concave (i.e. for (e pp , e sp ) = (e, e) and 
(1.33e, vL33e)) the diffusivity is approximately a constant. However, when the diffusivity 
shows exponential dependence on solvent concentration (for e pp = 2.0e) then the solvent 
concentration profile curve becomes convex. For the case in which the concentration is ap- 
proximately linear (e pp = e sp = 1.33e), the diffusivity increases linearly with concentration. 
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FIG. 4: [Color online] Solvent concentration profiles are plotted as a function of the scaling 
variable zt' 1 / 2 for (e pp ,e sp ) = (a) (e, e), (b) (1.33e, y/U&e), (c) (1.33e, 1.33e), (d) (2.0e,1.7e), and 
(e) (2.0e,2.0e). The red (light dark) solid lines represent the theoretical curve based on the solution 
of Eq. El 

Numerical solution of Eq. ^ for a given functional form of diffusivity on concentration has 
shown a similar relation between the shape of the concentration profile and the dependence 
of diffusivity on concentration.— 

B. Diffusion Coefficients from Equilibrated Polymer solution simulations 

Self- and corrected diffusion coefficients. The self-diffusion, D s (c), and corrected 
diffusion, D c (c), coefficients are calculated from Eq. El and Eq. [3 respectively. D s (c) and 
D c (c) as a function of solvent concentration for different polymer-polymer and solvent- 
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FIG. 5: Diffusivity D{c) as a function of solvent concentration calculated from the solvent concen- 
tration profile using Eq. The different symbols are for (e pp , e sp ) = (e, e) (closed circles), (1.33e, 
\/l-33e) (closed squares), (1.33e, 1.33e) (closed triangles), (2.0e, 1.55e) (open circles), (2.0e, 1.7e) 
(open squares), and (2.0e, 2.0e) (open triangles). The solid line in (b) is an exponential fit to the 
data. 

polymer interactions are shown in Fig. Efa) and (b), respectively. Note that the diffusion 
coefficients calculated for (e pp , e sp ) = (1.33e, vL33e) and (2.0e, v^Oe) is limited to low sol- 
vent concentration since the systems phase separate for large c. The critical value of solvent 
concentration above which the system phase separates can be approximately determined 
directly from the behavior of the solvent concentration profile shown in Fig. ^ and H3 We 
have observed that for (e pp , e sp ) = (1.33e, \/T33e), (2.0e, v^Oe), and (2.0e, 1.7e) the system 
phase separates for concentration values corresponding to the steep part of the concentra- 
tion profile. The critical solvent concentration value is approximately the concentration 
value corresponding to the inflection point of the concentration curve. The critical solvent 
concentration value for (e pp , e sp ) = (2.0e, y/2.0e) is basically the dilute limit. 

In general, the self and corrected diffusion coefficients shown in Fig. El show an expo- 
nential dependence on concentration. However, for (e pp , e sp ) = (e, e) and (1.33e, VL33e) the 
dependence of the corrected diffusion on concentration is weak at low solvent concentration. 
For a given value of concentration, as expected, both D s (c) and D c (c) decrease as the state 
of the polymer changes from melt to glassy. Note that D s (0) ~ D c (0) in all cases. 
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FIG. 6: Dependence of diffusion constants on solvent concentration, (a) D s {c) and (b) D c (c). 
Symbols are for (e pp , e sp ) = (e, e) (circles), (1.33e, \/1.33e) (squares), (2.0e, 1.7e) (triangles), and 
(2.0e, 2.0e) (stars). 

Diffusivity. To calculate the diffusivity D(c) using Eq. the thermodynamic factor 
given by Eq. |§]has to be first determined. Using the GCMD simulation method, the activity 
coefficient of the solvent is determined as a function of solvent concentration and is shown 
in Fig. 0on a ln-ln plot. Note that the concentration at which all the activity coefficients 
converge is the pure solvent case. As expected, the activity coefficient is constant for low 
solvent concentration and thus D(0) ~ D c (0) = D s (0) for all cases. As the solvent concen- 
tration increases the activity coefficient for (e pp , e sp ) = (e, e) and (1.33e, vL33e) decreases 
with solvent concentration and increases for e pp = e sp = 1.33e and 2.0e. This results in the 
diffusivity for the former two cases to be lower while for the latter two cases to be higher 
than the corrected diffusion constant. However, the curves in Fig. are not smooth, mak- 
ing it difficult to determine the thermodynamic factor and thus the diffusivity with high 
precession. 

Instead of numerically differentiating the activity curves as we did in paper I, the activity 
curve for each case was fit to a cubic function and the thermodynamic factor was determined 
analytically. The diffusivity calculated following this procedure is shown in Fig. |H1 For 
comparison, the diffusivity calculated from the solvent concentration profile reported in Sec. 
IVA is also shown as closed symbols. In general, the diffusivity calculated from the two 
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FIG. 7: Activity coefficient of solvent as a function of solvent concentration. Symbols are for (e pp , 
e sp ) = (e, e) (circles), (1.33e, vL33e) (squares), (2.0e, 2.0e) (triangles), and (1.33e, 1.33e) (stars) 

different approaches show good agreement. The diffusivity calculated using the Darken 
equation shows the expected behavior discussed above, the diffusivity is constant for e pp = e 
and increases exponentially for e pp = e sp = 2.0e. However, there is more scatter due to the 
uncertainty in the thermodynamic factor. 




o 0.2 

c 



FIG. 8: Diffusivity D(c) as a function of solvent concentration and open symbols are from Darken 
equation Eq. |1] and closed symbols are from the concentration profile using Eq. |21 Circles and 
triangles are for a pp — € sp — e, and 2.0e, respectively. 

In the interdiffusion study the behavior of D(c) is related to the shape of the solvent 
concentration profile. Similarly, in the present approach the behavior of D(c) can be directly 
related to the functional form of the activity coefficient. In general, for a constant diffusivity 
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the activity coefficient decreases with concentration, and increases with concentration for a 
diffusivity that increases with concentration. 

V. SUMMARY 

The effect of polymer-polymer and solvent-polymer interactions on the behavior of the 
interdiffusion of a solvent in to an entangled polymer matrix have been studied using large 
scale molecular dynamics and grand canonical Monte Carlo simulation techniques. By vary- 
ing the polymer-polymer interaction the state of the polymer is changed from melt to glassy. 
Correspondingly the solvent density profile changed from error function like to a sharp front, 
characteristic of Case II type transport, when Berthelot's rule is applied for the solvent- 
polymer interaction. The weight gain by the polymer matrix increased as t 1 ! 2 in agreement 
with Fickian diffusion, even for the case of a glassy polymer. This suggests that the precursor 
of the front is Fickian in agreement with recent experimental observations that characterize 
Case II diffusion by a sharp concentration front with Fickian type precursor.— i 15 i 18 i 19 The 
front, however, does not move in the time scale of our simulation. From simulation of 
equilibrated solvent-polymer solution it was found that the glassy system with Berthelot's 
rule applied for the cross term is immiscible except in the dilute limit suggesting that the 
front may not move in to the polymer matrix. Increasing the solvent-polymer interaction 
enhanced the solubility of the system without changing the nature of the diffusion process. 

The solvent concentration profiles have been fitted using the one-dimensional Fick's model 
of the diffusion process. The diffusivity, D(c), shows strong dependence on the state of the 
polymer. Far above the glass transition D(c) is approximately constant and then becomes 
concentration dependent as the polymer becomes glassy. The shape of the concentration 
profile and the behavior of D(c) is found to be directly related. The diffusivity is constant 
when the solvent concentration profile is concave, shows exponential dependence on solvent 
concentration when the solvent concentration profile is convex. 

The diffusivity as a function of solvent concentration was also determined using the 
Darken equation for simulations of equilibrated solvent-polymer solution. The diffusivity 
calculated using this approach is in good agreement with the diffusivity calculated from the 
solvent concentration profile. 

The advantage of the Darken approach for determining D(c) is that it requires much 



15 



smaller system sizes than for the direct simulation of the interdiffusion process. However, 
each solvent concentration c has to be determined separately and the simulation time re- 
quired to determine the D c (c) and the fugacity / are quite long compared to the interdiffusion 
studies. This is because the scaled concentration profiles superimpose even after relatively 
short times. 
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